Statistical Methods in Medical Research

The increase in methods focused on various types of survival outcomes has allowed practitioners to analyze data that are difficult or expensive to prospectively observe. Still, there are populations that are challenging to study. For example, obtaining a representative sample of couples attempting to become pregnant is difficult due to the dynamic nature of the population. This has led to an increase in the use of cross-sectional designs yielding backwards recurrent survival outcomes. In this paper, we consider the analysis of a survival outcome where subjects are observed if they are at-risk for a separate dependent survival outcome. The motivation for this problem is to determine which factors are associated with time-to-fertility-treatment (TTFT) among women currently attempting pregnancy in a cross-sectional sample. We propose appending a marginal accelerated failure time (AFT) model on TTFT with a conditional model on time-to-pregnancy (TTP) given TTFT to account for their dependence and avoid biases. We address challenges that arise due to the censoring of TTFT and the resulting increased computational complexity. The performance is validated via comprehensive simulation studies. We apply our approach to data from the National Survey of Family Growth to estimate the association insurance type has on TTFT, and estimate the impact of fertility treatment on TTP.

ObjectiveThe study of the potential intermediate effect of several variables on the association between an exposure and a time-to-event outcome is a question of interest in epidemiologic research. However, to our knowledge, no tools have been developed for the evaluation of multiple correlated mediators in a survival setting. MethodsIn this work, we extended the multimediate algorithm, which conducts mediation analysis in the context of multiple uncausally correlated mediators, to a time-to-event setting using the semiparametric additive hazards model. We theoretically demonstrated that, under certain assumptions, indirect, direct and total effects can be calculated using the counterfactual framework with collapsible survival models. We also adapted the algorithm to accommodate exposure-mediator interactions. Results and conclusionsUsing simulations, we demonstrated that our algorithm performs better than the product of coefficients method, even for uncorrelated mediators. The additive hazards model quantifies the effects as rate differences, which constitute a measure of impact, with applications that can be highly informative for public health. Our algorithm can be found in the R package multimediate, which is available in Github.

AO_SCPLOWBSTRACTC_SCPLOWRecent advances in sequencing technologies have allowed collection of massive genome-wide information that substantially enhances the diagnosis and prognosis of head and neck cancer. Identifying predictive markers for survival time is crucial for devising prognostic systems, and learning the underlying molecular driver of the cancer course. In this paper, we introduce -KIDS, a model-free feature screening procedure with false discovery rate (FDR) control for ultrahigh dimensional right-censored data, which is robust against unknown censoring mechanisms. Specifically, our two-stage procedure initially selects a set of important features with a dual screening mechanism using nonparametric reproducing-kernel-based ANOVA statistics, followed by identifying a refined set (of features) under directional FDR control through a unified knockoff procedure. The finite sample properties of our method, and its novelty (in light of existing alternatives) are evaluated via simulation studies. Furthermore, we illustrate our methodology via application to a motivating right-censored head and neck (HN) cancer survival data derived from The Cancer Genome Atlas, with further validation on a similar HN cancer data from the Gene Expression Omnibus database. The methodology can be implemented via the R package DSFDRC, available in GitHub.

Environmental exposures can regulate intermediate molecular phenotypes, such as gene expression, by different mechanisms and thereby lead to various health outcomes. It is of significant scientific interest to unravel the role of potentially high-dimensional intermediate phenotypes in the relationship between environmental exposure and traits. Mediation analysis is an important tool for investigating such relationships. However, it has mainly focused on low-dimensional settings, and there is a lack of a good measure of the total mediation effect. Here, we extend an R-squared (Rsq) effect size measure, originally proposed in the single-mediator setting, to the moderate- and high-dimensional mediator settings in the mixed model framework. Based on extensive simulations, we compare our measure and estimation procedure with several frequently used mediation measures, including product, proportion, and ratio measures. Our Rsq measure has small bias and variance under the correctly specified model. To mitigate potential bias induced by non-mediators, we examine two variable selection procedures, i.e., iterative sure independence screening and false discovery rate control, to exclude the non-mediators. We evaluate the consistency of the proposed estimation procedures and introduce a resampling-based confidence interval. By applying the proposed estimation procedure, we find that more than half of the aging-related variations in systolic blood pressure can be explained by gene expression profiles in the Framingham Heart Study.

The development of chronic disease is a long-term process that involves multiple endpoints, and few methods can assess the health benefits of a treatment regime over the disease course. Existing multi-state Cox models estimate survival risks by state over time, which are difficult to use when comparing the effectiveness of treatment regimes. A discrete-time split-state framework has been proposed, which divides disease states into substates by conditioning on past history. As this framework is both "memoryless" and "memorable", the time-specific transition parameters can be synthesized into summary measures, substate-specific life year (SSLY), multimorbidity-adjusted life year (MALY), and disease path. In this paper, based on this framework, we propose to investigate the causal effects of static and dynamic treatment regimes on health benefits over the entire disease course, under the assumptions of constant confounders from baseline and instantaneous effects of interventions on transition rates. Our method can identify the optimal treatment regime that generates the most benefits using MALY, and illustrate the mechanisms of treatment regimes affecting disease progression using SSLY and disease path. In the application, we evaluated the cardiovascular benefits of smoking cessation using data from the Atherosclerosis Risk in Communities (ARIC) study, where the course of heart disease was modeled in healthy (S0), at metabolic risk (S1), coronary heart disease (S2), heart failure (S3), and mortality states (S4). Compared to the regime "being a smoker in S0-S4", the MALY was 0.53 (95% CI: 0.21, 0.96), 6.10 (95% CI: 4.88, 7.19), and 4.34 (95% CI: 3.02, 5.47) years higher for the regimes "being a smoker in S0 and S1 and stop smoking if a person develops S2, S3, or S4", "no smoking in S0-S4", and "being a smoker at the start of intervention and stop smoking if age>65y", respectively. In summary, our method can evaluate the health benefits of treatment regimes over the disease course, and has the potential to improve the precision of chronic disease prevention.

ObjectivesIn clinical and epidemiological studies, the modified Poisson and least-squares regression analyses for binary outcomes have been used as standard multivariate analysis methods to provide risk ratio and risk difference estimates. However, their ordinary Wald-type confidence intervals can suffer from finite-sample biases in the robust variance estimators, and the coverage probabilities of true effect measures are substantially below the nominal level (usually 95%). To address this issue, new accurate inference methods are needed. MethodsWe propose two accurate inference methods based on the estimating equation theory for these regression models. A remarkable advantage of these regression models is that the correct models to be estimated are known, that is, conventional binomial regression models with log and identity links. Using this modeling information, we first derive the quasi-score statistics, whose robust variances are estimated using the correct model information, and then propose a confidence interval based on the regression coefficient test using{chi} 2 -approximation. To further improve the large sample approximation, we propose adapting a parametric bootstrap method to estimate the sample distribution of the quasi-score statistics using the correct model information. In addition, we developed an R package, rqlm (https://doi.org/10.32614/CRAN.package.rqlm), that can implement the new methods via simple commands. ResultsIn extensive simulation studies, the coverage probabilities of the two new methods clearly outperformed the ordinary Wald-type confidence interval when the regression function assumptions were correctly specified, especially in small and moderate sample settings. We also illustrated the proposed methods by applying them to an epidemiological study of epilepsy. The proposed methods provided wider confidence intervals, reflecting statistical uncertainty. ConclusionsThe current standard Wald-type confidence intervals may provide misleading evidence. Erroneous evidence can potentially influence clinical practice, public health, and policymaking. These possibly inaccurate results should be circumvented using effective statistical methods. These new inference methods would provide more accurate evidence for future medical studies.

In this article, we propose a Proportional Relative Hazards (PRH) model to differentiate subjects according to their risk for a primary event relative to competing events. The model estimates effects on the baseline ratio of the hazard for a primary event, or set of primary events, relative to the hazard for a competing event, or set of competing events ({omega}+ ratio). An analogous model is presented to estimate effects on the baseline ratio of the hazard for a primary event (or set of events) relative to the hazard for all events ({omega} ratio). A weighted regression method is introduced, along with practical presentation of risk-stratification using the PRH model in breast and head and neck cancer data sets.

Mediation analysis helps uncover how exposures impact outcomes through intermediate variables. Traditional mean-based total mediation effect measures may suffer from the cancellation of opposite component-wise effects, and existing methods often lack the power to capture weak effects in high-dimensional mediators. Additionally, most high-dimensional mediation analysis methods have focused on continuous outcomes, with limited attention to binary outcomes, particularly in case-control studies. To fill this gap, we propose an R2 total mediation effect measure within the liability framework that offers a clear and intuitive causal interpretation, provides additional insights beyond the mean-based measures, and is invariant to disease prevalence. We develop a cross-fitted, modified Haseman-Elston regression-based estimation procedure tailored for mediation analysis in case-control studies, which can also be applied to cohort studies. Our estimator remains consistent in the presence of non-mediators and weak effects, as demonstrated in extensive simulations. Theoretical justification for consistency is provided under mild conditions and without requiring exact mediator selection. In a case-control substudy of the Womens Health Initiative involving 2150 individuals, we found that many metabolites were mediators with weak effects in the path from BMI to coronary heart disease, and we estimated that 89% (95% CI: 57%-100%) of the BMI-explained variation in underlying CHD liability is mediated by the measured metabolomics. The proposed estimation procedure is implemented in the R package "r2MedCausal", available on GitHub.

Before applying flexible nonparametric models such as a generalized additive model (GAM), it is natural to ask whether a simpler parametric form suffices. To address this question, we develop TAPS (Test for Arbitrary Parametric Structure), a framework that integrates estimation and hypothesis testing to evaluate whether a prespecified parametric form adequately captures a covariate effect in a GAM. TAPS accommodates diverse structures, including linearity, piecewise linearity with changepoints, discontinuities with jumps, among many others. It is implemented in the R package mgcv.taps built directly on mgcv, enabling seamless adoption, broad outcome support, and scalability to biobank-scale data. Using UK Biobank data, we analyze 38 continuous and 8 binary traits to investigate two scientific questions: does the effect of a polygenic risk score (PRS) vary with age beyond a linear interaction, and does retirement at age 65 modify this age-varying effect? We find that age-varying PRS effects are common and often strongly non-linear, and that retirement at 65 significantly modifies these effects for five traits after multiple-testing correction.

AO_SCPLOWBSTRACTC_SCPLOWUnderstanding treatment effects on health-related outcomes using real-world data requires defining a causal parameter and imposing relevant identification assumptions to translate it into a statistical estimand. Semiparametric methods, like the targeted maximum likelihood estimator (TMLE), have been developed to construct asymptotically linear estimators of these parameters. To further establish the asymptotic efficiency of these estimators, two conditions must be met: 1) the relevant components of the data likelihood must fall within a Donsker class, and 2) the estimates of nuisance parameters must converge to their true values at a rate faster than n-1/4. The Highly Adaptive LASSO (HAL) satisfies these criteria by acting as an empirical risk minimizer within a class of cadlag functions with a bounded sectional variation norm, which is known to be Donsker. HAL achieves the desired rate of convergence, thereby guaranteeing the estimators asymptotic efficiency. The function class over which HAL minimizes its risk is flexible enough to capture realistic functions while maintaining the conditions for establishing efficiency. Additionally, HAL enables robust inference for non-pathwise differentiable parameters, such as the conditional average treatment effect (CATE) and causal dose-response curve, which are important in precision health. While these parameters are often considered in machine learning literature, these applications typically lack proper statistical inference. HAL addresses this gap by providing reliable statistical uncertainty quantification that is essential for informed decision-making in health research.

Mendelian randomization (MR) is a popular instrumental variable (IV) approach, in which genetic markers are used as IVs. In order to improve efficiency, multiple markers are routinely used in MR analyses, leading to concerns about bias due to possible violation of IV exclusion restriction of no direct effect of any IV on the out-come other than through the exposure in view. To address this concern, we introduce a new class of Multiply Robust MR (MR2) estimators that are guaranteed to remain consistent for the causal effect of interest provided that at least one genetic marker is a valid IV without necessarily knowing which IVs are invalid. We show that the proposed MR2 estimators are a special case of a more general class of estimators that remain consistent provided that a set of at least k{dagger} out of K candidate instrumental variables are valid, for k{dagger}[≤] K set by the analyst ex ante, without necessarily knowing which IVs are invalid. We provide formal semiparametric theory supporting our results, and characterize the semiparametric efficiency bound for the exposure causal effect which cannot be improved upon by any regular estimator with our favorable robustness property. We conduct extensive simulation studies and apply our methods to a large-scale analysis of UK Biobank data, demonstrating the superior empirical performance of MR2 compared to competing MR methods.

Compartmental models based on ordinary differential equations (ODEs) quantifying the interactions between susceptible, infectious, and recovered individuals within a population have played an important role in infectious disease modeling. The aim of the present paper is to explain the link between stochastic epidemic models based on the susceptible-infectious-recovered (SIR) model, and methods from survival analysis. We illustrate how standard software for survival analysis in the statistical language R can be used to estimate pivotal parameters in the stochastic SIR model in the very much idealized situation where the epidemic is completely observed. Extensions incorporating interventions, age structure and heterogeneity are explored and illustrated.

In real world data (RWD) studies, observed datasets are often subject to left truncation, which can bias estimates of survival parameters. Standard methods can only suitably account for left truncation when survival and entry time are independent. Therefore, in the dependent left truncation setting, it is important to quantify the magnitude and direction of estimator bias to determine whether an analysis provides valid results. We conduct simulation studies of common RWD analytic settings in order to determine when standard analysis provides reliable estimates, and to identify factors that contribute most to estimator bias. We also outline a procedure for conducting a simulation-based sensitivity analysis for an arbitrary dataset subject to dependent left truncation. Our simulation results show that when comparing a truncated real-world arm to a non-truncated arm, we observe the estimated hazard ratio biased upwards, providing conservative inference. The most important data-generating parameter contributing to bias is the proportion of left truncated patients, given any level of dependence between survival and entry time. For specific datasets and analyses that may differ from our example, we recommend applying our sensitivity analysis approach to determine how results would change given varying proportions of truncation.

A standard assumption for causal inference from observational data is that one has measured a sufficiently rich set of covariates to ensure that within covariate strata, subjects are exchangeable across observed treatment values. Skepticism about the exchangeability assumption in observational studies is often warranted because it hinges on investigators ability to accurately measure covariates capturing all potential sources of confounding. Realistically, confounding mechanisms can rarely if ever, be learned with certainty from measured covariates. One can therefore only ever hope that covariate measurements are at best proxies of true underlying confounding mechanisms operating in an observational study, thus invalidating causal claims made on basis of standard exchangeability conditions. Causal learning from proxies is a challenging inverse problem which has to date remained unresolved. In this paper, we introduce a formal potential outcome framework for proximal causal learning, which while explicitly acknowledging covariate measurements as imperfect proxies of confounding mechanisms, offers an opportunity to learn about causal effects in settings where exchangeability on the basis of measured covariates fails. Sufficient conditions for nonparametric identification are given, leading to the proximal g-formula and corresponding proximal g-computation algorithm for estimation. These may be viewed as generalizations of Robins foundational g-formula and g-computation algorithm, which account explicitly for bias due to unmeasured confounding. Both point treatment and time-varying treatment settings are considered, and an application of proximal g-computation of causal effects is given for illustration.

PURPOSEWe propose and demonstrate the feasibility and desirability of a novel model-free dose-finding design for phase I clinical trials. METHODSThe Generalized 3+3 (G3+3) design uses a set of simple rules summarized as follows: For 3 or 6 patients at a dose, apply the 3+3 design for making dosing decisions. For other numbers, if the observed toxicity rate (OTR) is less than 0.2, escalate to the next higher dose; if the OTR is greater than 0.29, de-escalate to the next lower dose; otherwise, stay at the current dose. RESULTSThe G3+3 design is the only design that can replicate the decisions of the 3+3 design for 3 or 6 patients among the popular designs compared like BOIN and i3+3. G3+3 generates desirable decisions when the number of patients treated is not 3 or 6, like the popular designs. Computer simulation verifies the superior operating characteristics of the G3+3 design. CONCLUSIONThe G3+3 design generalizes the popular 3+3 design so that desirable decisions can be made for any number of patients at a dose. G3+3 does not rely on statistical models, is simple and transparent, and can be implemented without a software tool. Therefore, it is expected to facilitate and enhance modern phase I dose-finding trials and early-phase drug development.

Adjustment for prognostic covariates increases the statistical power of randomized trials. The factors influencing increase of power are well-known for trials with continuous outcomes. Here, we study which factors influence power and sample size requirements in time-to-event trials. We consider both parametric simulations and simulations derived from the TCGA cohort of hepatocellular carcinoma (HCC) patients to assess how sample size requirements are reduced with covariate adjustment. Simulations demonstrate that the benefit of covariate adjustment increases with the prognostic performance of the adjustment covariate (C-index) and with the cumulative incidence of the event in the trial. For a covariate that has an intermediate prognostic performance (C-index=0.65), the reduction of sample size varies from 1.7% when cumulative incidence is of 10% to 26.5% when cumulative incidence is of 90%. Broadening eligibility criteria usually reduces statistical power while our simulations show that it can be maintained with adequate covariate adjustment. In a simulation of HCC trials, we find that the number of patients screened for eligibility can be divided by 2.7 when broadening eligibility criteria. Last, we find that the Cox-Snell [Formula] is a good approximation of the reduction in sample size requirements provided by covariate adjustment. This metric can be used in the design of time-to-event trials to determine sample size. Overall, more systematic adjustment for prognostic covariates leads to more efficient and inclusive clinical trials especially when cumulative incidence is large as in metastatic and advanced cancers. Key messagesO_LICovariate adjustment is a statistical technique that leverages prognostic scores within the statistical analysis of a trial. We study its benefits for time-to-event trials. C_LIO_LIPower gain achieved with covariate adjustment is determined by the prognostic performance of the covariate and by the cumulative incidence of events at the end of the follow-up period. C_LIO_LITrials in indications with large cumulative incidence such as metastatic cancers can benefit from covariate adjustment to improve their statistical power. C_LIO_LICovariate adjustment maintains statistical power in trials when eligibility criteria are broadened. C_LI

We extend Bayesian Logistic Regression to model the dose-toxicity relationship in the setting of phase I dose-escalation/ dose-finding trials for cancer immunotherapies. Immunotherapy drugs are associated with Cytokine Release Syn-drome, a systemic immune system reaction that can be mitigated when initial lower doses of the drug are administered to generate immune tolerance. This changes the classic dose-finding problem of determining an optimal safe dose, to a more complex problem where the search is for both the optimal safe dose and the dose regimen that allows patients to quickly and safely reach that dose without CRS. As part of solving this methodological challenge, we show how to jointly model CRS and non-CRS toxicities, which have distinct mechanisms, while controlling for the overall toxicity rate to make dose-escalation decisions.

There is a lack of approaches to evaluate the effectiveness of interventions when there are nonlinear impacts of covariates to the outcome series. Based on the classic framework of ITS/CITS segmented regression, while considering autocorrelation of time series, we adopted a nonlinear dynamic modeling strategy (Hammerstein) to measure the nonlinear effects of covariates, and proposed four optimized models: ITS-A, CITS-A, ITS-HA, and CITS-HA. To compare the accuracy and precision in estimating the long-term impact of an intervention between the optimized and classic segmented models, we constructed a sequence generator to simulate the outcome series with actual characteristics. The relative error with respect to the true value was the accuracy indicator, and the width of the 95% CI and the truth value coverage rate of the corresponding 95% CI are the precision indicator for model assessments. The relative error of impact evaluation in the four optimized models was 4.49 percentage points lower than that in the classic models, specifically ITS-A (14.34%) and ITS-HA (21.47%) relative to ITS (26.66%), CITS-A (16.57%), and CITS-HA (17.94%) relative to CITS (21.59%). The width of the 95% CI of point estimate of long-term impacts in the optimized models was 0.1261, which was expanded by 58.71% compared with 0.0875 for the classic model. However, the optimized models covered the true value in all test scenarios, whereas the coverage rates of the classic ITS and CITS models were 73.33% and 83.33%, respectively. The optimized models are useful tools as they can assess the long-term impact of interventions with additional considerations for the nonlinear effects of covariates and allow for modeling of time-series autocorrelation and lag of intervention effects.

Mediation analysis quantifies how an exposure affects an outcome through an intermediate variable. We extend mediation analysis to capture the cumulative effects of longitudinal predictors on longitudinal outcomes. Our proposed model examines how mediators transmit the effects of the current and previous exposure on the current outcome. We construct a least-squared estimator for cumulative indirect effect (CIE) and used three approaches (exact form, delta method, and bootstrap procedure) to estimate its standard error (SE). The estimator of CIE is unbiased with no unmeasured confounding and independent model errors between mediator model and outcome model at all time points, as shown in statistical inference and in simulations. While three SE estimates are numerically similar, bootstrap procedure is recommended due to its simplicity in implementation. We apply this method to Framingham Heart Study offspring cohort to assess if DNA methylation mediates the association of alcohol consumption with systolic blood pressure over two time points. We identify two CpGs (cg05130679 and cg05465916) as mediators and construct a composite DNA methylation score from 11 CpGs, which mediates for 39% of the cumulative effect. In conclusion, we propose an unbiased estimator for CIE. Future studies will investigate the missingness in mediators and outcomes.

BackgroundThe surge of treatments for COVID-19 in the ongoing pandemic presents an exemplar scenario with low prevalence of a given treatment and high outcome risk. Motivated by that, we conducted a simulation study for treatment effect estimation in such scenarios. We compared the performance of two methods for addressing confounding during the process of estimating treatment effects, namely disease risk scores (DRS) and propensity scores (PS) using different machine learning algorithms. MethodsMonte Carlo simulated data with 25 different scenarios of treatment prevalence, outcome risk, data complexity, and sample size were created. PS and DRS matching with 1: 1 ratio were applied with logistic regression with least absolute shrinkage and selection operator (LASSO) regularization, multilayer perceptron (MLP), and eXtreme Gradient Boosting (XgBoost). Estimation performance was evaluated using relative bias and corresponding confidence intervals. ResultsBias in treatment effect estimation increased with decreasing treatment prevalence regardless of matching method. DRS resulted in lower bias compared to PS when treatment prevalence was less than 10%, under strong confounding and nonlinear nonadditive data setting. However, DRS did not outperform PS under linear data setting and small sample size, even when the treatment prevalence was less than 10%. PS had a comparable or lower bias to DRS when treatment prevalence was common or high (10% - 50%). All three machine learning methods had similar performance, with LASSO and XgBoost yielding the lowest bias in some scenarios. Decreasing sample size or adding nonlinearity and non-additivity in data improved the performance of both PS and DRS. ConclusionsUnder strong confounding with large sample size DRS reduced bias compared to PS in scenarios with low treatment prevalence (less than 10%), whilst PS was preferable for the study of treatments with prevalence greater than 10%, regardless of the outcome prevalence. Key MessagesO_LIWhen handling nonlinear nonadditive data with strong confounding, DRS estimated by machine learning methods outperforms PS in scenarios with low treatment prevalence (less than 10%). C_LIO_LIHowever, if having linear data and small sample size data with strong confounding, we did not observe DRS outperformed PS even when treatment prevalence was less than 10%. C_LIO_LIOur results suggested that using PS performed better compared to DRS in tackling strong confounding problems with treatment prevalence greater than 10%. C_LIO_LISmall sample size increased bias for both DRS and PS methods, and it affected DRS more than PS. C_LI